!Même programme que ann-eqinttest2 mais avec recollement
!résolution de f(x) = 1- x*\int_0^1 dy f(x*y), petits x
!résolution de f(x) = 1- \int_0^x dy f(y), grands x
!-----------solution: f(x) = \exp(-x)
module fonctions
implicit none
integer, parameter :: points = 3000 !nombre de points d'évaluation des fonctions
integer, parameter :: maxit = 200 !nombre maximal d'itérations
double precision, parameter :: infty = 10.d0 !valeur numérique pour l'infini
double precision, parameter :: lambda = 0.9d0 !vitesse de convergence (=1.0) versus stabilité (=0.01...) (?)
contains

!------------------------------------------------------------------------
!Méthode de simpson pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine simptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
double precision :: wgt
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n==1) somme = h*somme/2.d0
wgt = 2.d0
if (n>1) then
	do i=1,n-1
	  if (wgt == 2.d0) then
	    wgt = 4.d0
	  else
	    wgt = 2.d0
	  end if
	  somme = somme + wgt*y(i)
	end do
	somme = h*somme/3.d0
end if
return
end subroutine simptab

!------------------------------------------------------------------------
!Méthode du trapèze pour calculer une intégrale connaissant uniquement
!quelques points discrets. Explication des variables:
!  y       : tableau de dimension (n+1) contenant les valeurs y(xi)
!  h       : largeur du pas d'intégration
!  n       : nombre de pas d'intégration
!  somme   : valeur de l'intégrale sur la base des points fournis
subroutine traptab(y,h,n,somme)
implicit none
double precision, intent(in), dimension(0:points) :: y
double precision, intent(in) :: h
integer, intent(in) :: n
double precision, intent(out) :: somme
integer :: i
somme = y(0)
somme = somme + y(n)
if (n==0) somme = 0.d0
if (n > 1) then
	do i=1,n-1
	  somme = somme + 2.d0*y(i)
	end do
end if
somme = h*somme/2.d0
return
end subroutine traptab

end module fonctions

!----------------------------------------------------------------------------------
!----------------------------------------------------------------------------------

program anneqinttest
use fonctions
implicit none
double precision, dimension(0:points) :: fout,fin,tmp,fout2
double precision :: h,err,err2
integer :: i,j,k,recol,div

h = infty/dble(points)
div = 5 !centrage du recollement
recol = int(0.1*real(points)) !largeur de la bande de recollement: +- recol autour de points/div
do i=0,points !condition initiale
	fin(i) = dexp(-(h*dble(i))**2)
end do
do j=1,maxit

!petits c: 1 - x*\int_0^1 dy f(x*y)
	do i=0,points/div+recol !valeurs de c
		do k=0,int(1.d0/h)
			tmp(k) = fin(nint(dble(i)*dble(k)*h))
		end do
		call traptab(tmp,h,int(1.d0/h),fout(i))
		fout(i) = 1.d0 - (dble(i)*h)*fout(i)	
	end do

!grands c: 1 - \int_0^x dy f(y)
	do i=points/div-recol,points
		call traptab(fin,h,i,fout2(i))
	end do
	fout2 = 1.d0 - fout2

!recollement
	fout(points/div-recol:points/div+recol) = (fout(points/div-recol:points/div+recol) + fout2(points/div-recol:points/div+recol))/2.d0
	fout(points/div+recol+1:points) = fout2(points/div+recol+1:points)
!--------

	err = sum(dabs(fin-fout))/dabs(sum(fout))
	fin = (1.d0 - lambda)*fin + lambda*fout
	err2 = 0.d0 	!erreur à la solution exponentielle
	do i=0,points
		err2 = err2 + dabs(fin(i) - dexp(-h*dble(i)))
	end do
	write(*,'(" Iteration: ",i5,", err.relative: ",e9.2,", |f_n(x) - exp(-x)|/points: ",e9.2)') j,err,err2/dble(points)
	open(unit=1, file='ann-eqinttest.txt')
	do i=0,points
		write(1,*) h*dble(i),fout(i)
	end do
	close(1)
end do

write(*,*) "Fin du programme."
read(*,*)
end program anneqinttest